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Abstract 

This paper further develops the zero-dimensional (0D) hydrodynamic coronal loop 
model “Enthalpy-based Thermal Evolution of Loops” (EBTEL) originally proposed 
by Klimchuk et al (2008), which studies the plasma response to evolving coronal 
heating. It has typically been applied to impulsive heating events. The basis of 
EBTEL is the modelling of mass exchange between the corona and transition region 
and chromosphere in response to heating variations, with the key parameter being the 
ratio of transition region to coronal radiation. We develop new models for this 
parameter that now include gravitational stratification and a physically motivated 
approach to radiative cooling. A number of examples are presented, including 
nanoflares in short and long loops, and a small flare. It is found that while the 
evolution of the loop temperature is rather insensitive to the details of the model, 
accurate tracking of the density requires the inclusion of our new features. In 
particular, we are able to now obtain highly over-dense loops in the late cooling phase 
and decreases to the coronal density arising due to stratification. The 0D results are 
compared to a ID hydro code (Hydrad). The agreement is acceptable, with the 
exception of the flare case where some versions of Hydrad can give significantly 
lower densities. This is attributed to the method used to model the chromosphere in a 
flare. EBTEL is suitable for general use as a tool for (a) quick-look results of loop 
evolution in response to a given heating function and (b) situations where the 
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modelling of hundreds or thousands of elemental loops is needed. A single run takes a 
few seconds on a contemporary laptop. 
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1. Introduction. 


Since the recognition in the 1970s that the magnetically confined solar corona is 
comprised of discrete loops, a great deal of effort has been devoted to modelling the 
temporal evolution of loop plasma. One approach is to solve numerically the one- 
dimensional hydrodynamic (ID hydro) equations of mass, momentum and energy 
conservation along a magnetic field line (or strand, or loop) in response to an imposed 
time-dependent heating function representing a flare or smaller heating event (e.g. 
Peres, 2000). Of importance is the ability of such models to generate “observables” 
that can be used to interpret coronal data (e.g. Bradshaw and Cargill, 2006; Bradshaw 
and Klimchuk, 2011). 

ID hydro models have two difficulties. One is the optically thick chromosphere at the 
lower boundaries. In principle this requires a full radiative-hydrodynamic treatment 
(e.g. McClymont and Canfield, 1983a, b) but one can attach a simple lower 
atmosphere that preserves the essential physics (e.g. Klimchuk et al., 1987; Antiochos 
et al., 1999). The second, and more significant, difficulty is the limitation imposed on 
the computational timestep by thermal conduction in the transition region (hereafter 
TR). In static equilibrium loops (e.g. Martens, 2010) the downward heat flux implies 
a temperature scale height ( Lj ) of under 1 km in the TR, and even shorter in hot 
flaring loops. Resolving this requires a fine grid, but when modelling thermal 
conduction the timestep scales as the smallest value of Lj, implying long run times. 

There is thus a need for simple and fast ways of modelling the coronal response to 
time-dependent heating. “Zero-dimensional” (0D) models, which average over the 
loop’s spatial dimension (Kuin and Martens, 1982, Fisher & Hawley, 1990, Kopp and 
Poletto, 1993, Cargill, 1994, Klimchuk et al., 2008, Aschwanden and Tsiklauri, 2009) 
accomplish this. In addition to providing “quick look” results, 0D models are useful if 
a loop is comprised of many hundreds or thousands of thin, thermally isolated, 
randomly heated strands (Cargill, 1994), which conventional ID hydro modelling still 
finds a large task. 

The success of ID and 0D models depends on handling correctly the exchange of 
matter between the corona, TR and chromosphere in response to a changing coronal 
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temperature. While the above OD models all address this to varying degrees (see 
Cargill et al., 2011, hereafter Paper 3), we base our discussion here on the work of 
Klimchuk et al., (2008: hereafter Paper 1). In Paper 1 we developed a OD model 
whose centrepiece was the explicit calculation of the enthalpy flux to and from the 
corona. The model, called EBTEL: “Enthalpy Based Thermal Evolution of Loops”, 
divides a loop into coronal and TR parts, the boundary being defined as where 
thermal conduction changes from a loss to a gain. Whether the enthalpy flux is into, 
or out of, the corona depends on whether the TR can radiate away the downward heat 
flux. If it cannot do this, then material is “evaporated” into the corona, whose density 
then increases (e.g. Antiochos and Sturrock, 1978). If the TR radiation cannot be 
powered by the downward heat flux, then there is a downward enthalpy flux, and the 
coronal density decreases (e.g. Cargill et al., 1995). The model was compared with 
several ID hydro simulations of an impulsively heated loop (starting each time with 
the same initial conditions), and gave reasonable agreement. 

EBTEL relies on three parameters, the most important of which is the ratio of the TR 
to coronal radiative losses. They govern both the initial equilibrium and how the loop 
cools after impulsive heating. It has become apparent through use of EBTEL, and 
attempts to benchmark EBTEL results against other known solutions of loop cooling, 
that the choice of this parameter in Paper 1 may only be appropriate for a long, 
tenuous loop. Thus, this paper seeks to put the EBTEL model on a firmer foundation 
and also looks at a broader range of loop evolution. The physical principles are 
unchanged, but we have undertaken a re-evaluation of the three key parameters 
(Section 2). The result is a model that now can follow with a good degree of accuracy, 
when compared with a ID hydro code, the evolution of loops over a wide range of 
lengths and temperatures (Section 3). In Paper 3 we will provide a comparison of OD 
models and sources of potential discrepancy with ID models. 

2. The models 

2.1 The governing equations of EBTEL 

The details of the model are discussed in Paper 1 and are restated briefly here. The ID 
energy equation is: 
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— = v(F + p)- dF ^ + Q- n 2 A(T) 
dt ds ds 


( 1 ) 


^ rj - T 

where v is the velocity, E = ^ , F = -k q T 5 ' 2 — is the heat flux, (9(f) is a heating 

f-1 05 

function that includes both steady and time-dependent components, A (T) = %T a is the 
radiative loss function in an optically thin plasma as defined in Paper 1, Equation (3), 
and .v is a spatial coordinate along the magnetic field. We have assumed that the flow 
is always subsonic and that gravity can be neglected from the viewpoint of the 
energetics. There is also an equation of state p = 2 nkT . 


For a corona loop of half-length L and a transition region of thickness / ( «L ), we 
define the boundary between corona and TR as the location where conduction 
changes from a loss to a gain (Vesecky et al., 1979). Integrating Eq (1) from the top of 
the TR to the top of the loop and enforcing symmetry boundary conditions, we find: 

~^ = -Vv 0 + F 0 + LQ - R c (2) 

y - 1 dt y-l 

where “overbar” denotes an averaged coronal quantity, subscript “0” denotes a 
quantity at the base of the corona (or top of the TR) and R c = n 2 A(T)L . Note that the 
heat flux and enthalpy flux play equivalent roles in providing energy to the TR. 


Integrating over the TR, and assuming the heat flux and flow vanish at its base, the 
pressure derivative and the heating can be eliminated since l « L, giving: 
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y-l 


P o v o + Aq +R tr - 0 


(3) 


where R, r is the integrated radiative TR losses. Eq (3) can then be combined with Eq 
(2) to give an equation for the coronal evolution: 


1 ‘S’-e-kfc+JU 


(4) 


y-l dt 

Note that conduction does not appear in Eq (4) emphasising its role as an energy 
redistribution mechanism as opposed to an energy loss. To calculate the density 
evolution, we adopt a similar approach to the mass equation, and in the corona find: 

dn _n qVq _ y-l 


dt 


L 2kT 0 Ly 


(F 0 + R tr ) . 


(5) 
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The average coronal temperature then follows from the equation of state: 

1 dT _ 1 dp 1 dn ^ 

T dt p dt n dt 

To solve the set of coronal equations (4) - (6) for the primary variables T , n and p , 
we need to define R tr , To and Fo in terms of the coronal quantities. The conductive 
losses are defined in terms of the loop apex temperature (T a ): 

Fq = — ( 2 / 7 ) A' () 7^7 ' 2 / L (Paper 1: Eq 20), so that there are three temperatures that 

characterise the corona: T , T a and To. T a and To are defined in terms of the primary 
variable T as C 2 =T / T a , C 3 =T 0 / T a . Finally, we define a third parameter: 
C, = R„./R c which modifies Eqs (4) and (5) to: 


=Q- R ;( l+ c,) 

y- \ dt L 

(7) 

dn _ H„v 0 __C 2 (r-l) ( } 

dt L C 3 2 kTLy 

(8) 


Specification of C/_.? is then required to solve Eq (6), (7) and (8). 

2.2 Modifications to the EBTEL model 

The original version of EBTEL described in Paper 1 used constant values of C 1 - 3 . 
However, as EBTEL has become more widely used, it is clear that further analysis of 
C 1-3 is required. Here we want to address three questions: (i) are the “baseline” values 
used in Paper 1 correct? (ii) how can EBTEL handle stratification due to gravity in 
long loops and (iii) how does EBTEL handle loop cooling when radiative losses 
strongly dominate? A secondary issue is how a general radiative loss function with a 
multiple power law can be included. We will show that these questions can be 
addressed by adopting a variable C/ within any given loop study, while implementing 
minor changes to C 2 and C,\ C/ is the most important parameter because the TR 
radiation regulates the value of the heat flux or enthalpy flux out of the corona that 
powers it, which in turn feeds into the coronal density and temperature. 

2.2.1 Assessment of parameters: equilibrium loops 


6 



Paper 1 used static equilibrium loop solutions to calculate C 1 . 3 . This is likely to be 
valid around the time of peak density when coronal radiative and conductive losses 
are comparable. Two approaches were considered. One had fixed values for all 
temperatures, namely Cj = 4, C 2 = 0.87 and C? = 0.5, and this was used to produce all 
the Figures in Paper 1. The second used a polynomial fit for Cj and C? over the 
temperature range 1-10 MK based on solutions of the equilibrium energy equation 
(Tables 1 and 2 of Paper 1). However the values of Cj and C3 quoted for short loops 
and T > 3 MK (Table 1 of Paper 1) are incorrect, and so the polynomial fit for C/ and 
C? present in early versions of the publicly available EBTEL code should not be used. 

We now reassess C/._? for static loops. First we use a simple power law radiative loss 
function and defer to Section (c) below how to include a multiple power law. We use 
the approximation A(T) =%T a = 1.95 10" 18 T~ 2B as our baseline power law, modified 
to A(T) =1.1 10" 31 r 2 below 10 5 K to avoid unrealistic losses at low temperatures. To 
calculate hydrostatic thermal equilibrium numerically, a Runge-Kutta method is used. 
T a and L are specified, and a double iteration calculates the base pressure and 
(constant) heating such that the appropriate boundary conditions are satisfied (T = T a 
and dT/ds = 0 at top of loop) for a given base temperature and vanishing base heat 
flux. This gives the familiar scaling laws between T a , L, Q and po. 

The derivation of C1.3 for general equilibrium loops relies on first calculating values 
for loops with no gravity and for a single power law form of A(T). The details of this 
are presented in Appendix A. There we first use the work of Martens (2010) to 
demonstrate analytically that Ci and C? are independent of all parameters except the 
slope of the radiative loss function for the case with no low temperature correction to 
A(T). Modifying A(T) at low temperatures shows that in the absence of gravity, C/.j 
may be taken as constants over a wide range of temperature and length (and hence 
pressure and heating rate). For L = 2.5, 5 and 7.5 10 9 cm and T a between 0.5-10 MK, 
C 2 and C? are roughly constant with values of 0.9 and 0.6 respectively. C/ varies a 
little more with Ta, but can be taken as approximately 2. We propose these as the 
“baseline” values of the constants and now discuss how they are modified to include 
other effects for equilibrium loops. 
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(a) Gravity: Here the main effect is that, while the TR radiation is driven by the 
downward heat flux, and so for a given coronal temperature is roughly fixed, the 
coronal radiation falls due to the stratification. Thus larger values of C/ = R„/R c can 
be expected for loops with significant values of the ratio A JL, where 
A = c 2 s lg = (2kTJ m )/ g is the gravitational scale height. We have thus solved the 

hydrostatic equations using the simple power law loss function for two loop lengths: 5 
x 10 9 and 7.5 x 10 9 cm and temperatures between 5 x 10 5 and 4 x 10 6 K. In the upper 
panels of Figure 1, the stars denote C; when gravity is absent (around 2 in all cases) 
and the circles C? when gravity is included. [Note that static solutions for T a = 5 x I 0 5 
K with L = 7.5 10 9 cm could not be found: see also Serio et al., 1981.] Cj increases as 
the temperature and scale height decrease. Note that for the loop shown in Figure 1 of 
Paper 1 (T a = 2-3 MK, L = 7.5 10 9 cm), Cj = 4 seem s reasonable. 


We now seek a parameterisation of the form Ci(T, L). There is little dependence of C/ 
on loop length itself, rather the key parameter is the ratio LA,. We write: 


c R,r T ^(g= Q ) T ^(g=°) ~ 

1 R c L* fr te = 0)l* c te = 0)i. R o 


(9) 


where quantities labelled “g=0” are values when gravity is neglected. The lower 
panels of Figure 1 shows the ratios R tr (g=0)/R tr (stars) and R c (g=0)/R c (circles) as a 
function of T a . The TR ratio is roughly unity which can be explained by noting that in 
static equilibrium, the energy balance in the TR is between the downward heat flux 
from the corona and the TR radiation [Eq (3)]. For given L and T a , the heat flux is 
roughly the same irrespective of whether the corona is stratified or not, so the TR 
radiation must also be roughly the same. The coronal ratio shows the expected drop in 
radiation when gravity is included. Now the second ratio in (9) is known to be 2, so if 
we have a simple expression for the third, we can calculate C; for different stratified 


loops. We argue that, for a given coronal temperature, 


R\g = 0) (n 2 T a ) g=0 (n 2 \ 


n 2 T a 


assuming the coronal length is the same with and without gravity and the averages are 
as defined in Paper 1, Eq (1 1). 
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Next assume that the coronal pressure is given by: p(s) ~ p 0 exp 


2 L . 

= sm 

xX 


7TS 

2 L 


for a 


semi-circular loop where we use a scale height based on the average temperature ( 
T = C 2 T a ) and that hydrostatic density stratification occurs only in the corona. We 
have calculated the average pressure by integrating p(s) numerically and find that this 
average value is well approximated by using the actual pressure at s/L = 0.4. So, the 
average pressure is written as: p = p 0 exp(-2Tsin(;r/ 5)/ Ax) to finally give: 


C, = R,r = R " (g 0) exp(4L sin(;r / 5) /Ax). (11) 

R c R c (g = 0) 

The plus signs in the upper panels of Fig 1 show that these values of Ci work 
reasonably well for all but the lowest temperature. 


(b) Multiple power law radiative loss function: Next, neglecting gravity, we 
evaluate Cy for a more complicated loss function by comparing results for the EBTEL 
loss function and the single power law one: 


K _ 

RJT) 

R tr (a = 

-2/3)" 

R c (cc - -2/3)" 


i 

m 

(N 

1 

II 

i 

_R c (a = 

-2/3) 

R e cn 


where R tr (T) and R C (T) evaluate the loss functions at a given temperature using the full 
power law in EBTEL. The right hand plot in Figure 2 shows little difference in the TR 
losses between the two radiative loss models (stars), so we can assume the first term 
in (12) is unity. The explanation is once again that the TR losses are determined by 
the heat flux from the corona. The coronal loss (circles) does show differences 
between the models. The second term in (12) is 2. For the third term, we can use the 
average coronal temperature (T -T = C 2 T a ) in Eq (12) to evaluate the radiation. The 

left hand plot of Figure 2 shows the same quantities as the lower panels of Figure 1 
for a loop of length 5 10 9 cm. This model for Cy is almost independent of the loop 
length. 

(c) Gravity and multiple power law loss function. We now combine the two 
corrections for loops with gravity and the general EBTEL loss function by replacing 
the ratio before the exponential in Eq (11) (which has no gravity and the simple loss 
function) with Eq (12) (which has no gravity and the full loss function), and using the 
fact that the TR losses are roughly the same for all cases: 
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— ^' ,a — =^— e xp(4Lsin(;r/5)/ A;r) 

JL Ug = 0,T) J 

(13) 

where we now denote the “equilibrium” value of Cj as C'/(eq). The first ratio will be 
taken as 2 in this paper, but is written in a general form to allow for changes to the 
coronal losses that no longer use our power law approximation. Figure 3 shows the 
results in the same format as Figure 1 for two temperature ranges and a loop length of 
5 10 9 cm. Other than at very low temperatures, (13) does acceptably. 

2.2.2 Radiative cooling phase 


_ R,(g,T) jR tr (g = 0,a = -2/3) 
R c (g,T ) \_R C (g -0,a--2/3) 


Equilibrium conditions are a good approximation if the strand (or loop) is evolving 
slowly, but may not be if it evolves rapidly. Consider a nanoflare event. The strand is 
far from equilibrium during the impulsive heating and the start of the thermal 
conduction cooling phase that follows. TR radiation is energetically insignificant at 
this time, so Cj does not affect the evolution. Static equilibrium is a reasonable 
approximation during the next phase of cooling when both thermal conduction and 
radiation, and possibly also enthalpy, are important. However, during the final phase 
of cooling when either radiation or enthalpy dominates, the equilibrium value of C; is 
a poor approximation. 


It is well known that for short, hot loops, there is a scaling T ~ n 2 ~ 2 5 during radiation- 
dominated cooling (Serio et al., 1991; Cargill et ah, 1995; Bradshaw and Cargill, 
2005, 20 10a, b), with a scaling approaching T ~ n for longer, more tenuous loops 
(Bradshaw and Cargill, 2010b). We can use the EBTEL equations adapted to radiative 
cooling to determine the appropriate value of C/. In the absence of gravity and 
neglecting thermal conduction and heating, Eq (4) and (5) are: 


1 dp _ 1 +R ) dn _ (7~ l )Rtr 
y-ldt L { c trh dt 2kT 0 Ly 

and so, on writing T ~ n 5 we can relate T and n as follows: 


(14) 


1 dp 1 dn _ T 0 (R c +R tr ) _ C 3 (1 + CQ 
p dt n dt T R tr C 2 Q 


( 15 ) 


This can be solved for Cj as: 
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-1 


(16) 


C j (rad) = 


C 2 (1+S) I 

^3 J 


where To is now the temperature at which enthalpy changes from a loss to a gain 
(Bradshaw and Cargill, 20 10a, b) and we denote C/(rad) as the value of C/ in the 
radiative phase. For 5 = 2(1), and the same values of C? and C? as for equilibrium 
loops, we find C/ = 0.59(1.25). Bradshaw and Cargill (2010b) suggest that 8 can be 
even larger for small loops, in the range 2.3 - 2.4 which gives C/(rad)= 0.5. We chose 
this as the baseline value for the radiative phase in the absence of gravity and 


a = -2/3. 


To include gravity and a full A(T), we adopt the same approach as in Section 2.2.1 on 
the basis of our work on radiative cooling (Cargill et al., 1995; Bradshaw and Cargill, 
2005, 2010a, b): 

C, (rad) - — - 0.5 ^ c ^ ex p(4^ s i n (7r / 5) / A tt) . (17) 

R c L R c (g = 0,T) 

Thus, including gravity increases C/ corresponding to the decrease in 8 discussed in 
Bradshaw and Cargill (2010b). 

2.2.3 Overall implementation of Ci 


We now wish to implement a formalism for C/ that takes on the equilibrium value 
when conductive and radiative losses are of the same order, and the radiative value 
when conduction is unimportant. We define the density for equilibrium conditions as 
n eq , and note that conductive (radiative) cooling of an impulsively heated loop at a 
given temperature is characterized by n < (>) n eq (e.g. Cargill, 1994, Cargill and 
Klimchuk, 2004). A formula of the generic form C/ = F(n/n eq ) is considered, and we 
require C; = Ci(eq) when n/n eq = 1 and C/ = C/(rad) when n/n eq is large. We define 
an “equilibrium” loop to be one where the coronal conductive loss is twice the coronal 
radiative loss (C/ = 2). The following form of F is considered: 


C, =F(n/n eq ) = 


C C Q+C,(rad)(nln eq ) m 

1 + (n/n eq r 


(18) 


where m is taken to be 1 or 2. C c o is the “conductive” value of C; and is adjusted to 
give Ci = Ci(eq) for n = n eq . One thus gets: 
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C c0 = 2 Q ( eq ) - C { (rad) 



In Paper 1 we calculated separate coronal and transition region differential emission 
measures (DEMs), the latter by two distinct methods. The DEM is defined as: 
DEM(T) = n 2 (dT /ds)~ l . Our modifications to EBTEL do not change the way the 


coronal DEM is calculated since the coronal parameters are our primary variables. On 
the other hand, the TR DEM relies on an assumption of constant pressure in the loop, 
which the introduction of gravity will invalidate. In Paper 1 we calculate the DEM by 
solving the following quadratic equation for dT /ds: 




A (T) = 0 


While Jo is constant and determined by the mass flow to and from the corona, the 
pressure in the last term is a TR quantity. Thus, when gravity is important we need to 
modify this term to account for the fact that the TR pressure will be larger than the 
coronal one. This is done by using our coronal pressure modification in reverse, so we 
write p TR = ]Iexp(2Tsin(;r/ 5)1 An) . This feature is demonstrated in the next Section. 
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In Paper 1 we also provided approximate form s of the DEM for three cases of loop 
evolution: strong conduction-driven evaporation, equilibrium, and strong radiative- 
driven condensation (draining). Of these, the third is unmodified, while the first two 
both involve the TR pressure, and need to be changed accordingly. 

2.3 The Hydrad Code 

An important part of our work is a comparison of the EBTEL results with the 1-D 
Hydrad code. This code has been fully documented in Bradshaw and Mason (2003a, b) 
so that we summarise the details briefly here. Hydrad solves separate time-dependent 
electron and ion energy equations together with equations of mass and momentum 
conservation and an equation of state. For the EBTEL comparison, we introduce an 
anomalously high electron-ion collision frequency that ensures 7) = T e under all 
conditions. For radiative losses, we supress the full capability of Hydrad to track 
multiple ion species, and use the same optically thin loss function as in EBTEL. 

In the figures, we show Hydrad results averaged over the top 80% of the loop. With 
the exception of the long loop (Case 3 below), there is little difference in the average 
T and n when the top 10, 50 or 80% of the loop is considered. For long loops the 80% 
average gives a density 25% larger than the 10% one, which can be attributed to the 
gravitational stratification present. 

3 Results 

We present four examples of loop heating. Each case is characterised by: (i) a loop 
half-length L, (ii) a background heating function which in turn implies a pre-event 
temperature and density, and (iii) a flare or nanoflare heating function with the form 
of a triangular pulse. Both the background heating and flare/nanoflare heating are 
uniform in space. 

3.1 Nanoflare heating of a short loop 
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The first case is a nanoflare in a short loop. A heating pulse with a half-width of 1 00 
sec and peak magnitude 10" 2 ergs cm" 3 s" 1 is released in a loop with half-length 25 
Mm. This has a total energy release of 2.5 10 9 A ergs cm" 2 , where A is the loop cross- 
sectional area. For a strand diameter of 200 km, we get a release of (rc/2)xl0 24 ergs. 
The background heating is 10" 5 ergs cm" 3 s" 1 giving an initial temperature and density 
of 0.53 MK and 1.3 10 s cm" 3 respectively. 

For clarity we show the new developments with EBTEL in two consecutive plots. The 
four panels of Figure 4 show the loop temperature and density as a function of time 
(top left and top right), the relationship between T and n (lower left) and C/ (lower 
right), for cases where neither gravitational nor radiative loss corrections are included 
in Ci. The plots show: (i) the EBTEL-1 results (C; = 4, C? = 0.89, Cj = 0.5: dots), (ii) 
the new “fixed” constant values (C; = 2, C;? = 0.9, C? = 0.6: dash-dot), (iii) the Flydrad 
simulation (thick solid), (iv) the linear C/ profile (m = 1 in Eq (18), dashed) and (v) 
the quadratic C/ profile (m = 2 in Eq (18), thin solid). The small panel in the 
temperature plot shows the evolution around the maximum. 

We start by noting that the temperature evolution for all the methods is broadly 
similar. C/ has little influence on the peak temperature since radiation is not playing 
any significant role at that time: the peak arises from a balance between heating and 
thermal conduction. EBTEL-1 has the lowest peak temperature (3.86 MK), followed 
by Flydrad (3.95 MK) and then EBTEL-2 (4.13 MK). The various models have 
slightly different initial densities for the same background heating which will have a 
small influence on the peak temperature. It can also be seen that after the decay phase 
there is a temperature undershoot for the quadratic C? model as it tries to recover the 
initial equilibrium. 

On the other hand, the density profiles do depend on the method used. EBTEL-1 has 
the lowest peak density of the EBTEL runs (1.1 10 9 cm" 3 ) and peaks soonest. If we 
compare this with 0 = 2 which peaks at 1.17 10 9 cm" 3 , the larger value of Ci requires 
the corona to lose more energy both by conduction and, more precisely, reduces the 
enthalpy transfer to the corona (the transition region radiation is stronger and can 
accommodate more of the heat flux). The two non-uniform models of Cj fall between 
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the first two cases. Hydrad gives a lower peak density of 10 9 cm" 3 and has oscillations 
due to material sloshing back and forward in the loop in response to the heating. 

The density in the decay phase shows more striking differences. Using a constant 
value of Ci = 2 increases the density and gives a (fortuitous, as we shall see) 
agreement with Hydrad. The linear and quadratic models increase the density further, 
with the discrepancy of the quadratic model being the more significant. The T-n 
scalings in the decay phase are also informative. To help the eye in this panel we 
include two thin solid lines showing T ~ n m (equilibrium) and T ~ n (radiative 
cooling) scalings. C; = 4 gives a T-n scaling corresponding to an “equilibrium” loop 
and examination of the losses from the corona due to conduction and radiation show 
that they are indeed equal after a few hundred seconds. A constant C; = 2 does not 
give the expected T-n scaling in the radiative phase. The non-uniform Q models give 
the T-n scalings that are closer to what is expected, with the quadratic profile giving 
the faster transition and a clearly over-dense loop. The Hydrad simulation gives a 
shallower slope in the T-n relation. 

However, just including a variable Cj to account for radiative cooling is not the whole 
story. In Figure 5 we explore the effect of including the corrections due to gravity and 
the multiple power law A(T). The figure has the same format as Figure 4, but we now 
show EBTEL-2 results using (i) the quadratic C/ (dotted), (ii) the correction to the 
quadratic case due to the full loss function (dash-dot), (iii) the correction due to 
gravity and the loss function for the quadratic profile (solid), linear profile (dashed) 
and the Hydrad results (thick solid). The radiative correction does not make a major 
difference, but the gravity correction does. The peak density moves slightly closer to 
the Hydrad value, but the agreement in the decay phase, where gravity becomes 
important, improves. The two models for Cj straddle the Hydrad solution, with the 
linear profile slightly low, and the quadratic one slightly high. 

The changes introduced by gravity can also be seen in the Ci plot. Coronal 
stratification is sufficiently important (L/A of order 0.4 or more after 1000 sec) that in 
the radiative phase one moves back towards values of C/ more typical of equilibrium 
loops rather than the small values typical of radiative cooling in flare loops. 
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In Paper 3 we discuss more fully the fact that getting a reasonable answer for the 
density in OD models is much harder than the temperature, the radiative phase in 
particular being awkward. We thus feel that the agreement in the density between 
EBTEL-2 and Hydrad shown here is satisfactory. This is only achieved by using a 
model for Ci that evolves in time and includes gravitational stratification. Finally, the 
fact that the linear and quadratic C/ models straddle the Hydrad result means we 
cannot yet comment on the preference of one over the other. 

Figure 6 shows the DEMs for some of these cases. These are computed over the 5000 
secs of the run, and so are averages over the life of the strand. On the left, the coronal 
DEM is shown for four runs: original EBTEL (dash-dot), quadratic C/ profde (dotted) 
and quadratic and linear C/ modified for gravity (solid and dashed). The thick dashed 
line shows the TR DEM for the gravitationally modified quadratic profile. In fact, for 
this case the gravitational correction to the TR DEM makes little difference. On the 
right are the total (corona + TR) DEMs for EBTEL- 1 (dash-dot), EBTEL-2 with 
corrections (solid) and Hydrad (thick solid). All the EBTEL models give reasonable 
agreement with Hydrad above 1.5 MK as would be expected from the previous 
figures. Below that temperature, the coronal value from EBTEL- 1 falls significantly 
compared to the EBTEL-2 results. However, when the TR part is included, the total 
DEM for both EBTEL models compares well with Hydrad: the only feature of note is 
a small depression in the EBTEL- 1 DEM compared the Hydrad and EBTEL-2 around 
1 MK. 

3.2 Nanoflare heating of a long loop 

The second case re-examines the long, tenuous loop discussed in Paper 1. The loop 
has a half-length of 7.5 10 9 cm and is heated by a nanoflare with amplitude 1.5 10° 
ergs cm" 3 s" 1 and half-width 250 secs. There is a background heating of 5 10~ 6 ergs cm" 
3 s" 1 giving an initial temperature and density of 0.52 MK and 3 10 7 cm" 3 respectively. 
When this example was presented in Paper 1, reasonable agreement with the 1-D 
hydro code ARGOS (Antiochos et ah, 1999) was noted. 
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Figure 7 shows results in the same format as Figure 4 for: (i) EBTEL-1 (dash-dot), (ii) 
EBTEL-2 with quadratic density profile (dots), (iii) as (ii) with gravity and radiation 
correction for quadratic and linear C; profile (solid and dashed) and (iv) Hydrad 
(thick solid). In (iii) the radiative correction has a small influence. As in the first 
example, the temperature throughout shows little difference between the four runs 
whereas the density varies considerably. We see that the unmodified quadratic C/ 
model gives very high densities in the decay phase that disagree significantly with the 
Hydrad and ARGOS results. However, modifying C/ to include gravity gives a 
density that is closer to the EBTEL-1 model. The Hydrad density is lower than all the 
EBTEL mns, the discrepancy being smaller at peak density than in Paper 1. In the 
decay phase, the linear Ci profile now does better than the quadratic one when the 
gravitational correction is included. Thus, we see that the agreement between EBTEL- 
1 and ARGOS in Paper 1 was partly coincidental. Interestingly the constant value of 
Ci = 4 used in EBTEL-1 does very well for this example. If one draws a straight line 
at Ci = 4 across the lower right panel, one could argue that this is indeed a reasonable 
first approximation to our variable C/. 

The DEMs are shown in Figure 8 in the same format as Figure 6. We see that 
inclusion of the quadratic Ci profile introduces overdense loops and so increases the 
DEM below the peak temperature but introduction of the gravitational correction 
removes this and we recover something quite similar to that in Paper 1. There is a 
small pressure correction to the transition region DEM as shown by comparing the 
thick dashed (includes effect) and dot-dashed (omits effect) in the left panel. 

3.3 Nanoflare with higher background heating 

Both the previous examples had a nanoflare with energy much larger than the 
background thermal energy in the loop or, equivalently, the background heating is 
much smaller than the peak nanoflare heating. However, nanoflare heating is not 
necessarily confined to a single heating / cooling cycle in a loop. Evidence now 
suggests that impulsive heating on occasions may be occurring in loops or strands that 
have not undergone such a cycle (e.g. Warren et al, 2011), so that the heating takes 
place in a higher ambient density. 
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We have re-run the case described in Section 3.1 with a higher background heating. 
To conduct a proper comparison with Hydrad, care is needed in setting up the initial 
state. The reader will have noticed in the earlier figures that the initial density and 
temperature was different in Hydrad and EBTEL for the same background heating. 
[This is a reflection of the fact that modelling the full spatial structure of a loop will 
give a slightly different value of the temperature and density from an approximate 
model that neglects this structure.] So long as the nanoflare is “large” in the sense that 
at its peak intensity the heating is much stronger than the background, this does not 
matter. However, for a “weak” nanoflare, it is important to ran EBTEL and Hydrad 
with roughly the same initial density. Thus, Hydrad is initialised with a background 
heating of 4.1 10" 4 ergs cm" 3 s" 1 giving a temperature and density of 1.36 MK and 8.31 
10 s cm” respectively. EBTEL has a heating of 3.93 10" 4 ergs cm" 3 s" 1 giving a 
temperature and density of 1.48 MK and 8.31 10 s cm" 3 respectively. To compare with 
EBTEL-1, we have to increase the background heating by 50% to ensure starting from 
roughly the same temperature and density. The nanoflare itself is as in Section 3.1. 

Figure 9 shows the evolution of the temperature and density for EBTEL-1, EBTEL-2 
with the gravitational and radiative corrections with both linear and quadratic C/ 
profiles and Hydrad. The agreement between the two approaches can be considered 
satisfactory with the linear C/ profile doing a better job with the density. Comparing 
with Figure 4, we see lower peak temperatures and, although the peak density is 
higher when the initial density is higher, the actual density increase is smaller. In the 
case of Section 3.1, the peak temperature is set by a balance between the maximum 
heating rate and thermal conduction losses. With the higher initial density we have 
here, the temperature never reaches this point. Instead, the peak temperature 
corresponds to an approximate equality between the total energy input of the 
nanoflare and the increase in thermal energy, p/(y-l), at roughly constant density. 

3.4 Small Flare 

Our fourth example is of a modest flare in a short loop with L = 25 Mm. The 
maximum heating is 2 ergs cm" 3 s" 1 , and the pulse half-width 100 secs, so the total 
heating per unit area is 5 1 0 1 1 ergs cm" 2 which for a loop diameter 20% of the half- 
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length gives a total energy input of 4 10 29 ergs. We neglect “thick target” heating 
which was discussed in Paper 1 and remains part of the EBTEL code. 

The four panels of Figure 10 show the loop temperature and density in the same 
format as above for the EBTEL- 1 values, the various profiled versions of Cj and two 
ffydrad runs. For this flaring case the gravitational correction to C/ is not important. 
Once again, the temperature plots for all the approaches are superficially similar, 
ffydrad gives the highest temperatures and EBTEL- 1 the lowest. The significant 
differences are seen in the density. The ffydrad density is much lower than all of the 
EBTEL runs: a peak of 2.5 10 10 from ffydrad contrasting with 4 10 10 from the various 
EBTEL runs, fn the decay phase the differences narrow, but we see from the T-n plot 
that while the scaling between T and n in the ffydrad run is consistent with our 
expectations, the EBTEL runs have a weaker decay (i.e. T n <2 ). 

We have addressed the density discrepancy from a number of viewpoints that are 
discussed more fully in Paper 3. These include: the coronal pressure forcing the 
chromosphere downward, leading in turn to longer loops and lower coronal densities 
(Klimchuk, 2006) and the violation of the subsonic flow assumption, ffowever, for 
this case, we believe that the low coronal density lies in the treatment of the 
chromospheric radiation in ffydrad. The baseline version of the code limits the 
radiative losses if the density above 30 kK exceeds 10 12 cm" 3 . This can be seen as 
modelling a number of processes, such as optically thick effects. The thick dashed 
line in Figure 10 shows what happens if this threshold density is set to 10 11 cm" 3 . Now 
there is much better agreement between EBTEL and ffydrad with the linear C i profile 
doing slightly better. The density limiter now does not permit the initial large 
downward heat flux to be radiated away at high density, but forces plasma to 
evaporate into the corona, as indeed EBTEL does. 

4. Discussion and Conclusions. 

Simple 0D hydrodynamic models have a long and quite successful history in 
modelling the temporal evolution of transiently-heated coronal loops. Motivated by 
discrepancies in some results, and by recent better understanding of coronal radiative 
cooling, we have updated our original version of the EBTEL model to include 
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gravitational stratification and correct radiative cooling. For the cases we have 
presented, we believe that the approximate model gives reasonable agreement with 
the full numerical results for the same parameters. EBTEL now permits highly over- 
dense loops to form in the cooling phase, and accounts for coronal stratification in 
long, cool loops. This is achieved through the evolution of the key parameter, the ratio 
of the transition region to coronal radiation. This is required to take on an equilibrium 
value around the time of maximum density, and a radiative value later on during the 
cooling. The transition is accomplished by either a linear or quadratic function. Our 
results suggest a slight preference for the linear one. 

EBTEL is a useful tool in looking at the generic evolution of temperature and density, 
as well as the DEM of single loops. It runs fast (a few seconds on a contemporary 
laptop), and can be convolved with other software to generate, for example, light 
curves in various coronal emission lines, DEM profiles as a function of temperature 
etc. But, perhaps more useful is the ability to model a multi-strand corona. In such a 
scenario (e.g. Cargill, 1994; Cargill and Klimchuk, 1997), the coronal emission comes 
from many (perhaps thousands) of separately evolving strands. This is still beyond the 
abilities of 1 D hydro codes, at least with a realistic tum-around time whereas EBTEL 
can model such a scenario in a few hours, and indeed perhaps less if a properly 
optimised version is used. 
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Appendix. Comparison of analytical and numerical values of Cj for a simple 
power law radiative loss function 


Ci and C 3 can be calculated analytically from Martens (2010). Assuming uniform 
heating, a single power law slope of a for the radiation function, and boundary 
conditions of vanishing heat flux at top and bottom of loop, and vanishing 
temperature at bottom, he writes the energy equation in terms of the variable 
P = (T/Tj 12 as: 


d 2 r/ 

ds 2 


£ = 


= ri M ~Z, M = - 


2(2 -a) 
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where T a is the apex temperature and the scaling laws have been used to eliminate L, 
Q and p. B(a,b) is a beta function. Martens further solves the energy equation for a 
variable u = p~ M as: 


s/L = 0 r (u, A + 1,1/2), 

where |3 r is the normalised incomplete beta function. 


At the point where conduction changes from a gain to a loss, denoted by subscript 


zero, (Al) gives r/ = 


3 + 2 a 


2(2-a) 


J 


or, in real units, — : 
T 


3 + 2a 


( 2 -a) 


Setting a =-1/2, we get C 3 = To/T a = (2/7 ) 2/5 = 0.606. For o= -2/3, C 3 = 0.584 


We can also calculate Ci as follows. The dimensionless coronal radiative losses are: 
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Now Eq (Al) integrates once to give, on application of the boundary conditions: 
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The total radiative loss is just £ in these units so that the TR loss is then: 


Silo) 


J 7l u ds = ^(Tj 0 ) + e 


f drj^ 

V ds y s=s(%) 


We have calculated r|o above, and so can obtain C'/, which is independent of Q, L and 
p. For a= -1/2, we get C; = 1.76 and for a= -2/3, Ci = 2.095. 


We now compare the Martens solutions with a numerical solution with a lower 
boundary at 2 x 10 4 K and a single power law loss function above 10 5 K and a loss 
function scaling as T 2 below. [This eliminates the problem that the vanishing heat flux 
is only exactly enforceable in the limit of vanishing base temperature.] Care is needed 
with the grid on which the hydrostatic equations are solved. We used: 

si L = (2/^)[sin' 1 x- Wl-x 2 

where x is evenly distributed between 0 and 1. 5000 points are used. The motivation 
for this grid can be seen from Eq (Cl) of Rosner et al (1978) and it does give well- 
resolved solutions at all temperatures. 

An array of cases has been run: three loop half-lengths, 2.5, 5 and 7.5 x 10 9 , and T a 
between 10 6 and 10 7 for each length. It turns out that the results are by and large 
independent of the loop half-length, so individual cases are not shown, rather the 
range of values found as shown in the following Table. 


a = -1/2. 

C 1 

c 2 

c 3 

Analytic 

1.76 

0.89 

0.606 

Numerical 

1.65-1.74 

0.895 

0.62-0.61 

a = -2/3 




Analytic 

2.09 

0.89 

0.585 

Numerical 

1.88-2.06 

0.892 

0.61-0.59 


Table Al. The constants Cj, C 2 and Cj for two loss functions. The range of values in 
each box are those obtained as T a increases from low to high. 
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It can be seen that C 2 = 0.9 and C 3 = 0.6 are reasonable values for both cases. The 
lower values of C? correspond to smaller T a where the T 2 loss function at lower 
temperatures makes a greater relative contribution to the loop losses. We would argue 
that for a simple model, Ci = 1.7 for a = -1/2 and C/ = 2 for a = -2/3 are appropriate. 
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Figure 1: The upper two panels show C? as a function of T a for L = 5 10 9 cm (left) 
and 7.5 10 9 cm (right) for a single power law loss function with a low temperature 
correction. Stars, circles and plus signs are, respectively, Ci in absence of gravity, C/ 
with gravity (both are from numerical solutions of the hydrostatic equations) and the 
estimate of Ci in Eq (1 1). The lower two panels show the ratio of the radiative losses 
without gravity to those with gravity in the transition region (stars) and corona 
(circles). The ratio of the two transition region losses is roughly constant. 

Figure 2: Ci as a function of T a for L = 5 10 9 cm. There is no gravity. In the left plot 
stars, circles and plus signs are, respectively, C? for single power loss function with 
low temperature correction, C/ for the full EBTEL loss function, and the estimate of 
Ci in Eq (12). The right column shows the ratio of radiative losses assuming a single 
power law and the full EBTEL form in the transition region (stars) and corona 
(circles). 


Figure 3: The upper row shows Ci as a function of T a for two temperature ranges and 
L = 5 1 0 9 cm. Stars, circles and plus signs are, respectively, C/ for single power loss 
function and no gravity, for the EBTEL loss function and gravity, and the estimate of 
Ci in Eq (13). The lower row shows the ratio of radiative losses assuming a single 
power law with low temperature correction and no gravity, and the EBTEL loss 
function and gravity in the transition region (stars) and corona (circles). 


Figure 4. The temperature and density (upper row) and Ci (lower right panel) as a 
function of time for a nanoflare in a short ( L = 25 Mm) loop. The lower left panel 
shows the relationship between T and n, and the evolution can be tracked in time by 
going right along the horizontal lines at low density (the heating phase) and following 
the curves. The line coding is as follows: EBTEL-1 values of C 1.3 (dot), EBTEL-2 
values of C1.3 implemented as follows: EBTEL-2 with constant C; (dash-dot), linear 
Ci model (dashed) and quadratic Cj model (thin solid). The Hydrad results are shown 
by the thick solid lines. The small panel in the temperature plot shows its evolution in 
more detail around the maximum. In the T-n plot, the straight solid lines are the 
scalings T~ n and T~n I/2 . Ci is only shown for the linear and quadratic models and 
the gravitational and radiative corrections to C? are not included. 
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Figure 5. As Figure 4, but shows the effect of introducing radiative and gravity 
corrections to the C/ model. The dotted curves show results for the quadratic C/ 
profile in Figure 4, the dash-dot curve includes only the radiative correction and both 
radiative and gravitational corrections are shown for linear (dashed) and quadratic 
(solid) cases. The thick, solid curve is the Hydrad result. 


Figure 6. The DEM for a nanoflare in a short (L = 25 Mm) loop. The left panel shows 
separate coronal and TR contributions. The dash-dot, dotted, solid and dashed lines 
are respectively the coronal DEM for (i) EBTEL-1 values of C/_j, (ii) the quadratic C/ 
profile without gravity and radiation corrections, and (iii) and (iv) quadratic and linear 
Ci profiles with gravitational and radiative correction included. The thick dashed line 
is the transition region DEM associated with coronal DEM (iii). The right panel 
shows the sum of the coronal and TR DEMs for EBTEL-1 (dash-dot), EBTEL-2 with 
quadratic transition and radiative and gravity modifications (solid) and Flydrad (thick 
solid). 

Figure 7. The results for a nanoflare in a long loop ( L = 75 Mm). The figure format is 
the same as Figure 4. The curves are: EBTEL-1 (dash-dot), EBTEL-2 with (a) 
quadratic C/ model (dots), (b) quadratic model with gravity and radiative correction 
(solid), (c) linear model with corrections (dashed) and Flydrad (thick solid). 

Figure 8. The DEM for a nanoflare in a long (75 Mm) loop. The line coding is as in 
Figure 6 except that the thick dash-dot line in the left panel is the TR DEM in the 
absence of the gravitational correction. 


Figure 9. A nanoflare in a small loop with large initial density. The dash-dot, solid, 
dashed and thick solid curves show EBTEL-1, quadratic and linear Cj profiles for 
EBTEL-2 and Flydrad. 

Figure 10 Temperature, density, T-n and C? as a function of time for a small flare. 
The line coding is as follows: dash-dot (EBTEL-1), quadratic Q transition, no 
corrections (dots), quadratic and linear C; models with gravity and radiation (solid 
and dashed). Two Hydrad results are shown corresponding to the two chromospheric 
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radiation models in the text. The one where the radiation is limited above a density of 
10 12 cm" 3 is shown as the thick solid line in all panels, that where the radiation is 
limited above a density of 10 11 cm" 3 is shown as the thick dashed line in the density 
plot. 
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